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An overset grid thin-layer Navter- Stokes code has been extended to include dynamic motion of helicopter 
rotor blades through relative grid motion. The unsteady flowfieid and airloads on an AH-1G rotor in forward 
flight were computed to verify the methodology and to demonstrate the method’s potential usefulness towards 
comprehensive helicopter codes. In addition, the method uses the blade’s first harmonics measured in the flight 
test to prescribe the blade motion. The solution was impulsively started and became periodic in less than three 
rotor revolutions. Detailed unsteady numerical flow visualization techniques were applied to the entire unsteady 
data set of five rotor revolutions and exhibited flowfieid features such as blade vortex interaction and wake roll- 
up. The unsteady blade loads and surface pressures compare well against those from flight measurements. 
Details of the method, a discussion of the resulting predicted flowfieid, and requirements for future work are 
presented. Overall, given the proper blade dynamics, this method can compute the unsteady flowfieid of a 
general helicopter rotor in forward flight. 


Introduction 

T HE accurate computation of a helicopter flowfieid is es- 
sential for proper and efficient airload predictions in for- 
ward flight and hover. The constantly changing aerodynamic 
environment and loads are important features of rotorcraft 
aerodynamics. Strong tip vortices in the rotor wakes dominate 
the flowfieid to produce a highly unsteady and nonuniform 
induced velocity field at the rotor disk. Recent emphasis on 
numerical simulation procedures for complex flows have made 
it possible to accurately obtain rotor blade aerodynamic char- 
acteristics by solving the governing differential equations. 

Descriptions of various computational fluid dynamics (CFD) 
methods for rotorcraft problems based on the solution of full- 
potential, Euler, and Navier-Stokes equations can be found 
in the literature. Notable among them is the introduction of 
embedded or overset grid scheme to the existing single grid 
scheme by Duque and Srinivasan. 1 Also, Srinivasan et al. 2 * 3 
used a thin-layer Navier-Stokes method for forward-flight 
simulation of a nonlifting rotor blade as well as rotor in hover. 
Chen et al. 4 and Agarwal and Deese 5 solved the Euler equa- 
tions. 

Unlike fixed wing aircraft, the helicopter airloads depend 
greatly upon the unsteady dynamic motion of the blades. To 
model the blade motion, comprehensive helicopter analysis 
methods have been developed that use aerodynamic and wake 
models coupled to structural dynamic models. The aerody- 
namics are typically based upon lifting-line theory and depend 
upon either linear methods or two-dimensional experimental 
airfoil data with corrections for unsteady and three-dimen- 
sional effects. 
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Notable among the comprehensive helicopter analysis codes 
is CAMRAD/JA by Johnson. 6 To account for compressibility, 
three-dimensional effects, and arbitrary blade geometries, 
Strawn and Tung 7 coupled CAMRAD to the full potential 
rotor (FPR) code. Strawn et al. 8 demonstrated the capability 
to model rotors in forward flight using the coupled CAMRAD 
and FPR by computing the airloads of the Puma rotor. 

Hernandez and Johnson 9 used the coupled CAMRAD/JA- 
FPR method to compute the AH-1G rotor in forward flight. 
They were able to compare well against flight test data. In 
addition, they found that tip core size specifications have a 
large effect upon the predicted rotor loads. 

Ramachandran et al. 10 developed a method based upon the 
vorticity-embedding technique. This method solves the un- 
steady full potential equation on a Eulerian grid with an over- 
set Lagrangian vortical velocity field. He also applied his method 
to the forward flight of the AH-1G rotor system with favorable 
comparisons to flight tests. 

All of the above methods use single-structured grids. Single 
grids have limited use. For example, it would be very difficult 
to use a single-structured grid to investigate the aerodynamic 
interference between either a rotor and fuselage or the aero- 
dynamic interference between a rotor and the rotor hub. In 
addition, almost all of the above methods use wake models 
to include the influence of the wake. Wake models limit the 
method’s generality and require considerable adjustments to 
their approximations to obtain accurate solutions. 

The work summarized herein uses an alternative griding 
method known as overset grids or the Chimera method 11 to 
allow for blade motions and to accurately compute the vortical 
wake from first principles. The overset-grid scheme greatly 
simplifies arbitrary blade motions, which is important in 
achieving trimmed flight conditions, as well as in predicting 
airloads accurately in an effort to develop a comprehensive 
helicopter analysis method. Also, with overset grids one can 
more easily discretize the domain with simple well-defined 
grids that accurately compute the rotor wake. 

Until recently, few methods were able to solve the transient 
flow about multiple bodies moving relative to one another. 
Meakin 12 computed the unsteady flowfieid of the V-22 tiltro- 
tor aircraft, including the rotor rotation, by solving the un- 
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Table 2 Computed and flight test rotor force 
and moment coefficients 



C 7 , 


C*f v , 

Cmv> 


thrust 

torque 

longitudinal 

lateral 

Computed 

0.00455 

0.000227 

-0.00022 

0.000008 

Flight 

0.00464 

0.0002196 

0.0 

0.0 


ft = 0.19. M r = 0.65, Re - 9.73 x 10*, and cr = 0.0651. 



Fig. 11 Time history of rotor thrust. 

revolution at various radial locations. The figures show strong 
blade vortex interaction for the advancing blade around azi- 
muthal locations of 70-90 deg and on the retreating side 
around 270 deg. The loads are underpredicted for the ad- 
vancing side. At the retreating side, loads are slightly under- 
predicted for the cases presented here. Blade-vortex inter- 
actions at outboard stations seem stronger as indicated by the 
negative lift. 

Table 2 lists the computed rotor force and moment coef- 
ficients. As shown, the rotor is not trimmed in pitch, but 
essentially trimmed in roll. The power is overpredicted by 
approximately 15%. 

Figure 11 illustrates the time variation of the integrated 
rotor thrust for the fifth rotor revolution. The solid line pre- 
sents the predicted azimuthal variation while the dashed line 
represents the time averaged or gross rotor thrust. The gross 
thrust agrees within 1.8% of the gross weight of the aircraft 
in flight. 

Conclusions and Future Work 

The unsteady forward-flight flowfield of the AH-1G heli- 
copter’s two-bladed rotor system was computed by solving 
the thin-layer Navier-Stokes equations on moving overset 
grids. The resulting flowfield visualizations and airload com- 
parisons verify the method’s versatility for rotorcraft prob- 
lems. The method efficiently obtained periodic solutions within 
three rotor revolutions. Unsteady streak lines showed the 
significant blade vortex interactions and wake roll-up. In ad- 
dition, the streak-line patterns showed that grid resolution 
needs improvement in the rotor wake to reduce diffusion in 
the wake. Blade surface pressures compared well with the 
flight test data, especially at the retreating side. This com- 
parison was somehow surprisingly better in the inboard sec- 
tion. Leading-edge surface pressure irregularities point to er- 
rors in the original airfoil coordinate definitions. The spanwise 
and integrated load data also correlates fairly well with the 
flight airload data and previous computations. 

Although the rotor was only partially trimmed, these results 
demonstrate the method’s capabilities. To trim the rotor prop- 
erly, one would need to include blade dynamics by tightly 
coupling the unsteady load prediction to a suitable dynamics 
model, which is a straightforward extension to the present 
formulation. This coupling would also allow for aeroelastic 
effects such as torsion and bending. The inclusion of a fuselage 
through additional overset grids would further improve the 
solutions. 
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Fig. 9 Computed and measured blade surface pressure coefficient 

distributions, rlR = 0.97; A. flight; , computation, = a) 90, 

b) 105, c) 270, and d) 300 deg. 

blade roots descend further than those from the tips. This 
behavior illustrates the wake roll-up process. It also gives 
more insight into the wake diffusion as the particles continue 
to spread in the far-field background grid. This region has no 
intergrid boundary points; therefore, the wake diffusion is 
inherent to the solution and can only be controlled by ade- 
quate grid resolution. 

The fuselage, hub, and other control surfaces have a sig- 
nificant effect on the blade loads. Although this simulation 
neglects those important components, the overall load and 
surface pressure comparisons have fair agreement with the 
flight test. The following load comparisons give considerable 
confidence in the computations. 

Calculated and measured pressure coefficients are pre- 
sented in Figs. 8 and 9 for different azimuthal and radial 
locations. Figure 8 shows the pressure distribution at the in- 
board radial station 0.6. The solid line is the computed so- 
lution while the symbols represent the experimental value. 
Note that the scales are different for each pressure plot. To- 
wards the upper surface leading edge, the inboard stations’ 
pressure distributions either overshoot or underpredict the 
pressure, while the i// = 180-deg position compares fairly well 
with a difference in the trailing-edge pressures. 


Figure 9 gives the surface pressure distributions at a radial 
station of rlR - 0.97. Again, the solid line is the computed 
solution while the symbols represent flight data. The lower 
surface pressures are well predicted for both advancing and 
retreating sides, except at 300 deg. The upper surface ad- 
vancing and retreating sides have a consistent trend. For the 
advancing blade, the upper surface consistently underpredicts 
the leading-edge pressure towards the leading edge. On the 
retreating side, the upper surface pressures compare well with 
the experiment with some small overshoots. Hernandez and 
Johnson 9 found similar behavior in their predicted pressure 
distributions. 

At all radial stations, the calculated pressure distributions 
depict irregular behavior towards the leading-edge region; 
especially at 90 deg in Figs. 8 and 9. A preliminary analysis 
of the airfoil coordinates of the associated surface slopes and 
inviscid two-dimensional airfoil calculations suggests that er- 
rors in the input blade geometry’s leading-edge region causes 
these irregularities. The same effects appear in the results of 
Hernandez and Johnson, 9 who used the same input data. This 
discrepancy needs to be corrected in future work. 

Figure 10 compares the unsteady measured (dashed line) 
and calculated (solid line) airloads of the fifth computed rotor 






.Q.2 1 i — 

0 90 180 270 360 

d) Azimuth 

Fig. 10 Computed and measured blade normal force coefficients 

showing the blade-vortex interaction; — , flight; , computation. 

r/R = a) 0.60, b) 0.75, c) 0.97, and d) 0.99. 
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Fig. 7 Computed flow visualization using unsteady streak-line particles. Particles are released at inboard and outboard blade tip every S deg 
of azimuthal rotation. Traces shown for total five rotor revolutions. Flow periodicity is approximately established after second revolution: a) top 
view and b) perspective view. 

moment and out of balance in pitching (longitudinal) hub 
moment. The important point is that the overall method is 
stable, it predicts airloads that agree well with the flight test, 
and it lays the foundation for more complete trim computa- 
tions by including an appropriate dynamic model. 

Figure 6 shows one of the blade’s normal force coefficient 
time history at r/R ~ 0.97. The solid line is the computation, 
whereas the dashed line is the experiment. The figure shows 
that the normal force has an initial transient that exists for 
the first rotor revolutions and then becomes periodic. 

Figures 7a and 7b illustrate unsteady particle streak-line 
patterns after five rotor revolutions. The unsteady flowfield 
analysis tool (UFAT) by Lane 22 at NASA Ames Research 
Center computed these unsteady streak lines. Rakes of par- 
ticles were released every' 5 deg of rotation along a plane of 
points just behind each blade’s tip and root trailing edge . Each 
rake has a different shading to distinguish them from each 
other. 

Figure 7a illustrates a top view of the highly complex flow- 
field caused by the wakes interacting with the blades. The 
particles released from a leading-blade cross over the follow- 
ing blade, showing evidence of blade vortex interaction. The 
particles also show evidence of excessive wake diffusion. Ide- 
ally, the tip vortices should maintain their core radius and 
intensity. The computed streakline particles remain close to 
each other initially, but when they encounter the following 
blade they tend to diffuse. This behavior indicates that the 
overset grid interpolations may be affecting the solution. 

With the present grids, once the wake leaves the blade it 
must pass through five intergrid boundaries (blade grid to 
intermediate grid to background grid back to intermediate 
and then another intermediate to blade). The interpolation 



Fig. 8 Surface pressure coefficient rfR = 0.60; ▲, flight; , com- 

putation. 

process tends to degrade the solution because of the required 
interpolation at the boundaries and the disparate grid reso- 
lutions. The number of intergrid boundary interpolations can 
be reduced by either changing the intermediate grid positions 
or consolidating the two separate intermediate grids into one 
with more effectively distributed points. 

The streak-line patterns in Fig. 7 also demonstrate the rotor 
wake descent and roll-up. The particles released from the 
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Fig. 4 Blade grid with intergrid hole boundary. 



Fig. 5 Global background grid with blade system. 


To improve interpolation, each blade gnd lies within an 
intermediate grid. The intermediate grids are Cartesian with 
points concentrated at the blade’s vicinity as shown in Fig. 4. 
The grids extend from the hub rotation axis to approximately 
five blade chords from the upper surface, seven chords from 
the lower, four chords in front, seven chords behind, and five 
chords from the tip. Each grid contains 75 points spanwise, 
75 points chord, and 41 points from bottom to top. Figure 4 
illustrates boundary planes within an intermediate grid with 
the blade grid cutting through it. 

The global background grid shown in Fig. 5 completes the 
overset grid system. The global grid extends to four blade 
radii from the hub center upstream, downstream, and to the 
sides. The grid also extends two blade radii above the blade 
and two and one-half radii below. The grid consists of a total 
of 95 x 95 x 51 points with points clustered vertically in the 
rotor disk vicinity. 

The entire moving overset system totals 1.62 million grid 
points. During the grid motions, the background grid remains 
stationary as the blade and intermediate grids rotate together 
through it. Figure 5 highlights the relative positions of the 
intermediate and blade grids with respect to the background 
grid. The intermediate grids move only in rotation about the 
spin axis and subsequently create hole regions within the back- 
ground grid. The blade grids also rotate about the spin axis, 
but then also pitch and flap about those respective axes. Dur- 


ing their pitch and flap motions the blade grids create holes 
within the intermediate grid as shown in Figs. 3 and 4. 

Flight Test Conditions 

The flight conditions chosen as a validation case comes from 
the flight tests documented by Cross and Watts 20 and Cross 
and Tu. 21 This test was chosen because of its extensive load 
survey, acoustic measurements, and detailed blade harmonics 
data. In addition, select flight conditions have been computed 
by two other numerical methods. 9 10 Both methods showed 
good correlations with the flight test data. 

The flight tests were performed with the AH-1G helicopter 
at NASA Ames Research Center. The rotor is a two-bladed 
rectangular-planform teetering rotor with the operational loads 
survey (OLS) symmetrical airfoil section, and a linear twist 
of - 10 deg from root to tip. Each blade has an AR of 9.8. 

The forward-flight case chosen is test point no. 2157 re- 
ported in Cross and Watts. 20 This flight condition has an ad- 
vance ratio of 0.19, hover tip Mach number of 0.65, flight 
Reynolds number of 9.73 x 10 8 , and a rotor thrust coefficient 
equal to 0.00464. These conditions correspond to a forward 
speed of 82 kn at a rotor rotation rate of 315.9 rpm. 

The first blade harmonics reported by the flight test were 
used to prescribe the rotor blade motions. Table 1 lists the 
blade harmonics. 

Results and Discussion 

The results were computed on the Numerical Aerodynamic 
Simulation Facility’s Cray C-90 computer located at NASA 
Ames Research Center. The time-accurate calculation im- 
pulsively starts from freestream conditions with the viscous 
no-slip boundary condition applied at the blade surfaces. The 
time steps for this rotor blade system correspond to approx- 
imately 0.3125 deg of rotation or a nondimensional time step 
of 0.082. Each rotor revolution requires about 1152 time steps. 
The complete unsteady computation required a total of 45 h 
of single processor CPU time for five complete rotor revo- 
lutions producing approximately 40 Gbytes of flowfield data. 

As noted in the previous section, the computed results used 
the first blade harmonics as measured by the flight test. How- 
ever, the blade collective was corrected to the value recom- 
mended by Ramachandran et al., 10 which allowed them to 
match the measured value of thrust. No further trim of the 
rotor was attempted in the present investigation with the result 
that the computed thrust was about 1 .8% too low. In addition, 
the rotor was approximately balanced in rolling (lateral) hub 


Table 1 Blade first harmonics, 
test point 2157 

ft.. e u , e u , £ lr , 

deg deg deg deg deg 

64) -5.5 L7 -0.15 Z13 

M = 0.19. M t = 0.65, Re = 9.73 x 10 ft . and 
C r = 0.00464. 



Fig. 6 Time history of normal force coefficient. 
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Fig. 2 Moving overset grid schematic. 


In overset grids, the quality of the solution interpolation at 
the boundaries depends on the relative grid cell ARs, skew- 
ness, and clustering, etc. It also depends on the boundary’s 
proximity to high flow gradient regions. While forming the 
grids, one must ensure that the boundaries have adequate 
overlap, nearly equal cell ARs, and low skewness. Boundaries 
should also avoid high flow gradients. 

With moving overset grids, individual grids move with their 
appropriate grid motion. As the grids move, the holes and 
hole boundaries change with time. Figure 2 provides a sche- 
matic for a helicopter blade in rotation, with one blade grid 
omitted for clarity. As the blades rotate from some time state 
T 0 to another time T n , the grid attached to the blades rotate 
along with them. Subsequently, the holes change with the 
blade rotation as shown. 

To determine the grid’s changing connectivity and hole points, 
the code known as domain connectivity functions in three 
dimensions (DCF3D) by Meakin 12 was employed. DCF3D 
uses inverse mapping of the computational space to limit the 
search time and to compute hole and outer-boundary inter- 
polation stencils. The major expense in DCF3D is the creation 
of the inverse maps. However, the maps are independent of 
the relative orientation of the grids and so it repeatedly uses 
the maps during the grid movement. 

During the flowfield solution process, intergrid boundaries 
are constantly created due to the grid movement. After each 
flow solution time step, grid connectivity data has to be de- 
termined. Currently, the flow solver (Chi-TURNS) and the 
connectivity algorithm (DCF3D) are separate Fortran codes. 
A UNIX C-Shell script is used to join the two methods. 

With each time step, DCF3D obtains the latest connectivity 
data and hole points. Relevant intergrid boundary data are 
written to disk. The hole-fringe points are then checked to 
ensure that they are introduced to the solution gradually. The 
flow solver then begins the solution process by reading in the 
DCF3D files from disk, solving the equations of motion, and 
then moving the grids to the next time step. Once complete, 
the solution and grids are written to disk and the whole process 
continues until the desired number of time steps have been 
completed. 


Blade Motion 

The method assumes rigid blade motions in flap and pitch. 
The periodic blade motion for pitch and flap as a function of 
blade azimuth can be described by a Fourier series 10 17 as 
shown in Eqs, (7) and (8): 

Pitch 'Vv- 

6 - 6 0 + 0, r cos tp + 0 ls sin ip + cos 2ip 

4- 0^. sin 2s + • - * (7) 

Flap 

P = ft, + p ]c cos </r -I- /3 b sin ip 4- cos 2i P 

+ sin 2s 4- •• * (8) 

Using only the mean and first blade harmonics, Eulerian 
angles prescribe the blade motion to the flow solver. Euler 
parameters or Eulerian angles 18 are useful and convenient 
ways to express motion of rotating bodies in terms of the fixed 
inertial frame. 

In this method, the blade rotates about its spin axis at some 
given rotational rate. At each time step, the blade rotates 
through by an increment of t p that results in a change in pitch 
and flap. The incremental change in the blade position is then 
imposed by transforming the position vector through succes- 
sive matrix multiplications as shown in Eq. (9): 

T = [A][B][C] x m = Tx 0]a (9) 

The transformation matrix T consists of the rotation ma- 
trices A, 5, and C. The matrices A, £, and C represent the 
various coordinate rotations. See Amirouche 18 for details of 
the transformation matrices. 

Grid System 

The grid system consists of five overset grids; one for each 
rotor blade and one intermediate transition grid for each blade 
to help connect them to the single global background grid. 
The rotor grid is of C-H topology with clustering near the tip, 
root, and leading and trailing edges. Beyond the tip and in- 
board of the root chord the surface grid collapses to a slit and 
forms a singular plane that extends approximately one chord 
length beyond either side. Each rotor blade grid has 119 total 
chordwise points and 45 points in the spanwise direction that 
lie on the body. In the normal direction, the first grid point 
is 0.00003 chords off the body and extends approximately 0.75 
chord from the surface in all directions. The hyperbolic grid 
generator by Chan et al. 19 generated the resulting volume grid 
illustrated in Fig. 3. 



Fig. 3 Blade grid cutting through intermediate grid. 
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steady thin-layer Navier-Stokes equations on moving overset 
grids. To accomplish this calculation, he developed a method 
that determines the interpolation coefficients between the 
overset grids. His domain connectivity algorithm (DCF) was 
verified in detail in his work. But the set of flight conditions 
he considered were hypothetical and did not include the dy- 
namic pitch, flap, or lag motions of the rotor blades. The 
method presented here accounts for the blade motion by ap- 
propriately moving the overset grid system according to the 
blade harmonics measured in a flight test. 

In summary, a numerical method that computes the un- 
steady flowfield of a helicopter rotor in forward-flight or hover 
is presented. The method solves the thin-layer Navier-Stokes 
equations on a system of moving overset grids and uses the 
blade harmonics measured in flight to prescribe the blade 
motions. Comparisons between computed and flight-test un- 
steady blade loads and surface pressures of the AH-1G heli- 
copter are presented. Numerical flow visualizations exhibit 
flowfield features such as blade-vortex interaction and wake 
roll-up. This article provides details of the method, a discus- 
sion of the resulting predicted flowfield, and demonstrates 
the capability to model the unsteady flowfield of an arbitrary 
helicopter rotor in forward flight. 

Methodology 

Algorithm 

The single-grid flow solver developed by Srinivasan et al. 2 
with overset-grid modifications implemented by Duque and 
Srinivasan 1 was employed in the study. Calculations by Sri- 
nivasan and Ahmad 13 have demonstrated the method’s ver- 
satility by computing the hovering flowfield of a rotor and 
whirl stand. The method solves the thin-layer Navier-Stokes 
equations shown in Eq. (1): 


Background Grid 



Fig. 1 Schematic of overset grid system. 


The matrices L and V are formed by performing either 
backward or forward differences on the appropriate flux Ja- 
cobians A , B, or C as shown in Eqs. (3). In Eqs. (3), At is 
the time step, A f and represent forward and backward 
differences, respectively, I is the identity matrix, and a is the 
spectral radius. The matrix D is a diagonal matrix that com- 
pletes the back-solve process. This decomposition scheme re- 
duces the required floating point operations in comparison to 
block tridiagonal methods: 

L = I + IBj k ,\t(-A /Jc , + V^A„ — B j k I 4 - ^ v B n 
~ Cm, + 

D = / + lB hkJ M(A- + B n + C- + A; 

+ b: + c;) hkJ U) 

V = / + IB^jbtiAfrj + A f A- + Bfa 
+ & V B~ + C+ kJ + A^C~) 


-Q + -E + —F + —G = Re' 1 — S (1) 

dt di dri d{ d( ' ' 


where Q = J~ ! (p, pp, pv , poj, e) is the vector of the conserved 
variables, density, momentum, and energy, scaled by the 
transformation Jacobian, and £, F, G, and 5 are the scaled 
in viscid and viscous flux vectors. 

With overset grids a sequence of grids are placed such that 
they lie arbitrarily within a primary grid. For example, Fig. 
1 provides an example where an airfoil curvilinear grid lies 
within a background Cartesian grid. The airfoil grid captures 
features such as the boundary layers, tip vortices, and shocks, 
etc. The background grid surrounds the airfoil grid and carries 
the solution to the far field. The background grid was gen- 
erated with some knowledge of the airfoil’s surface and outer 
boundary locations. Consequently, some of the background 
grid points lie within the airfoil's solid body regions and must 
be removed from the solution. Once removed, hole regions 
remain within the interior of the larger background grid and 
create a set of boundary points known as hole-fnnge points. 
The airfoil grid interpolates data to the background grid at 
the background’s hole-fringe points. Conversely, the back- 
ground grid interpolates data to the airfoil grid at the airfoil 
outer boundary points. 

Within each separate overset grid, different flow solver 
methods may be used to solve the governing equations of 
motion. The current method uses an implicit upwind method 
for all of the individual grids. The method advances the so- 
lution in time using the lower-upper symmetric Gauss-Siedel 
(LU-SGS) implicit scheme by Yoon and Jameson 14 shown 
in Eqs. (2). The scheme is third-order accurate in space and 
first-order accurate in time: 


c . Q* = - IBAtL ~ X( 3i 

Solving 

LD 'UAQ = -IBAtVi i Q** = DQ* (2) 

AQ = U-'Q ** 


where A - = A ± <j x . 

With overset grids, the additional array IB marks the hole 
regions within the grid interior by taking on the value of 0. 
Outside of the hole and in the valid regions of the flowfield 
the IB array equals 1. The IB array removes the hole regions 
from the solution set as shown in Eqs. (2). 

The flux terms use a Roe upwind-biased scheme for all 
three coordinate directions with higher-order MUSCL-type 
limiting to model the shocks accurately. 15 The resulting flux 
differences are shown in Eq. (4): 

& = V ( E* + A ( E~ + V n F; + A „F- + V.G; + A ( G~ 

(4) 

At the hole fringe points the flux evaluations reduce to 
second order at the resulting interior boundaries by modifying 
the primitive variable evaluations with the IB array. Equation 
(5) shows the ^-coordinate flux evaluations: 

E(Q L , Q k ) - «£{&,) + E(Q l ) - \A(Q l , Q„)\(Q k - QJ] 

(5) 

The primitive variable quantities Q are obtained bv first- 
order extrapolation functions at boundaries or third order 
away from boundaries. As shown in Eqs. (6), these quantities 
drop to first order through the use of the IB array: 

Ql = Q, + IHQ, - £?>-,) + «e>., - 
Qh = Q,+ , - [*((?/., - Q, ) + KG,*, - G>. 

( 6 ) 

Finally, the flow solver assumes fully turbulent flow for the 
blades and in viscid flow in background grids. The simple al- 
gebraic turbulence model of Baldwin and Lomax 16 is used to 
estimate the eddy viscosity. 



